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We propose a model for the fluctuation dynamics of the local denaturation zones 
(bubbles) in double-stranded DNA. In our formulation, the DNA strand is model as 
a one dimensional Rouse chain confined at both the ends. The bubble is formed when 
the transverse displacement of the chain attains a critical value. This simple model 
effectively reproduces the autocorrelation function for the tagged base pair in the 
DNA strand as measured in the seminal single molecule experiment by Altan-Bonnet 
et. al (Phys. Rev. Lett. 90, 138101 (2003)). Although our model is mathematically 
similar to the one proposed by Chatterjee et al. (J. Chem. Phys. 127, 155104 (2007)) 
it goes beyond a single reaction coordinate description by incorporating the chain 
dynamics through a confined Rouse chain and thus considers the collective nature 
of the dynamics. Our model also shows that the autocorrelation function is very 
sensitive to the relaxation times of the normal modes of the chain, which is obvious 
since the fluctuation dynamics of the bubble has the contribution from the different 
normal modes of the chain. 



I. INTRODUCTION 

In 1953 Watson and Crick [1] proposed the structure of DNA to be a stable double- 
stranded helix. The stability comes through the staking interaction and the hydrogen bond- 
ing between the base pairs in the opposite strands. But actually this picture represents 
the equilibrium structure of DNA under physiological conditions. As the interaction energy 
between these base pairs is only few fc^T, even at room temperature due to thermal fluc- 
tuations locally DNA strand opens up creating what is called "bubbles". These bubbles 
have different sizes and lifetimes. The creation and annihilation kinetics of these bubbles is 
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termed as the breathing dynamics. On increasing the temperature or changing the pH these 
bubbles add up to form larger bubbles and eventually the double stranded DNA denatures. 

In recent past single molecule fluorescence correlation spectroscopy (FCS) experiment 
by Altan-Bonnet et. al Q has provided the first quantitative insight to the relaxation 
kinetics of the breathing mode of the double-stranded DNA. This breathing mode refers 
to local denaturation and reclosing of the double-stranded structure. In their experiment, 
two bases of the double-stranded DNA, corresponding to opposite strands are tagged with a 
fluorophore and a quencher respectively. So when the DNA structure is closed, fluorophore 
and the quencher are in close proximity and the fluorescence is quenched. But due to 
thermal fluctuation when the DNA structure opens up creating a bubble, the fluorescence is 
restored. Hence the base pair fluctuation leads to a fluctuation in fluorescence intensity. This 
fluctuation in fluorescence intensity is monitored by FCS which determines the characteristic 
dynamics of the relaxation of dynamic correlations in the fluctuation of base pairs. They 
introduced a correlation function 

G(t) = mm - (m (im m 

(i(oy) - (i(o)f 

where I(t) is the fluorescent intensity at time t. G(t) was found to be multiexponential. 
Interestingly the correlation function for all the DNA constructs, at all temperatures follow 
the same universal temporal behavior and when presented as a function of rescaled time they 
all collapse into a single universal curve G(t/ti/ 2 ), where t\/ 2 is such that G(t/ti/ 2 ) = 0.5. 
To explain the experimental data they proposed a simple kinetic model in which G(t/ti/ 2 ) 
becomes 
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G ex (z) = (1 + — )erfc( s j— ) - sj—e~- (2) 

where, z = and m is a parameter which is adjusted to 0.328 to ensure that G ex (l) = 
0.5. 

Since this experimental work by Altan-Bonnet et. al, the topic "bubble dynamics" as 
it is commonly refereed to has received a great deal of theoretical attention. Just few 
years after this seminal experimental study of the transient time-dependent rupture and 
re-healing of double stranded DNA, Bicout and Kats [9] proposed a kinetic scheme based 
on two state model (closed or open) of the double stranded DNA. Their formulation gave 
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an analytical expression for the survival probability, correlation function and life time for 
the bubble relaxation dynamics but formulation did not consider the structure of the double 
stranded DNA at any level. Very recently Srivastava and Singh proposed a theoretical 
model where the interaction between the base pairs of the opposite DNA strands is described 
by Peyard-Bishop-Dauxois (PBD) [12| potential and the separation between the base pairs 
( "y" in their notation) follows a Fokker-Planck equation. Thus only making the interaction 
realistic but still not taking into account of the dynamics of the chain and restricting to 
a single relevant dynamical variable (separation "y") description. In a paper by Jeon, 
Sung and Ree [10|] double stranded DNA was modeled as a duplex of semiflexible chains 
mutually bonded by weak interactions in other words they used an extended worm-like chain 
model and a Langevin dynamics simulation was performed to examine the size distribution 
and dynamics of the bubble. Shortly after this Fogedby and Metzler jsj came up with a 
theoretical model for the bubble dynamics . They used the Poland- Scheraga free energy 4] 
where the free energy is a function of bubble size x. The dynamics of x follows a Langevin 
equation and the corresponding Fokker-Planck equation is analogous to the imaginary time 
Schrodinger equation for a particle moving in a Coulomb potential subject to a centrifugal 
barrier. This mapping enabled them to calculate the correlation function. But the best fit to 
the long time behavior was obtained only when they considered the dynamics of x in a linear 
potential. Thus to reproduce the long time data one can merely start with a linear potential 
but unfortunately it may not fit to the experimental data in the intermediate time range. 
Interestingly the shortcomings of the Fogedly-Metzler model was pointed out by Chatterjee 
et. al j^j]. They assumed that the distance between fluorophore and the quencher follows 
an overdamped Langevin equation in a harmonic potential. Although their theory predicts 
the experimental data reasonably well it does not take into account of the fluctuation of 
different modes of the DNA strand contributing to the bubble dynamics. A better theory 
should account for those and actually our model does. In our formulation the DNA-strand 

fin 

is described by a Rouse chain [6|, [7| confined at both the ends. Transverse displacement of 
the string accounts for the bubble formation. Naturally our model takes into account of the 
contribution of different modes of the string to the breathing dynamics. 
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II. OUR MODEL 

We describe the bubble by a confined Rouse chain in a harmonic potential, V(x) = 
^K>R(n) 2 , where R(n,t) denotes the position of the nth segment of the confined Rouse chain 
in space and t is the time. In other words the chain is described by a field R(n, t). Because of 
the thermal fluctuations the bubble undergoes Brownian motion and its time development 
is described by the equation 

^dRin.t) ,d 2 R(n,t) „. , „, . 

(^Qf 1 = k ^ } - *R(n,t) + f(n,t). (3) 

In the above, ( is the friction coefficient for the nth segment and k is the force constant 
for the confining potential which accounts for the staking interaction between the opposite 
strands of the DNA molecule. 

The fluorescent intensity at time t, denoted by I(t) can be written as 

I(t) = A0(R(a,t) -a) (4) 

where R(a, t) denotes the position of the ath segment in space. Later on we will choose a 
such that it corresponds to the center of the bubble. As mentioned earlier in the experiment 
one measures the following correlation function. It is worth mentioning that this is the the 
relevant dynamical coordinate/variable in our model. Obviously it is not a one dimensional 
phenomenological reaction coordinate/dynamical variable like the separation between the 
donor and the acceptor as considered by Chatterjee et al. 5j but actually a collective 
dynamical variable in the sense that it has the contribution from all the normal modes of 
the chain as is shown later in Eq.(6). 

G(t) = (mm - (m (no)} 

So in order to evaluate G(t) we should first evaluate (I(t)I(0)) and (I(t)). In our model 
the correlation function, (I(t)) would be 



(1(f)) = A(6(R(a,t) -a)) 

Then we use the fact that the step function can be expressed as an integral over delta 
function and in the next step we use the fourier integral representation of the delta function. 
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Now R(n, t) follows Rouse dynamics and it has a boundary conditions R(0, t) = and 
R(L, t) = 0, where L is the chain length. Keeping these boundary conditions in mind one 
can express R(n, t) in terms of fourier modes 

oo 

R(n,t) = 2j2X p (t)Mn) (6) 
P =i 

4>(n) should satisfy the boundary conditions 0(0) = and 4>(L) = 0. X p is the pth normal 
mode of the chain. 

Substituting R(n,t) from Eq. (jSJ) into the expression for (I(t)) to get 

oo 
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The above average is evaluated as follows 
First we introduce 

b(s) = 2k 1 <j) p {a)8(t - s). 
Then the quantity in the angular bracket can be written as 

oo oo oo 

/ 2ik l Y,X p (t)Mn)\ / *E J" dsX p (s)fc(s)\ 



= 1 — OO 



Then we use the definition of the characteristic functional |8| and write the above quantity 
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as 



oo oo oo 



'E / dsX p (s)b(s)\ -1 / dti / dt 2 E 6(ti)H _1 (*i-*2)6(*2) 

p— 1 —oo \ — — p — oo — oo p— 1 



where 



^ 1 (ti-t 2 ) = (X p (t 1 )X p (t 2 )) 
With this the integral in the exponent becomes 



CXJ OO m 

oo 

dti / dh^bitjH-^h-hWh) = 4fc?C(t) 

-oo — oo ^ 
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where 

oo oo 

C(0) = £ P (a) 2 (X p (t) 2 ) = £ p (a) 2 (* p (0) J 
p=i p=i 

where we have used the fact (X p (t) 2 ) = (X p (0) 2 ) . Then 

oo oo 

and subsequently (I(t)) can be evaluated by integrating over a-y and k±. 

.00 
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Now the above expression for (/(£)) does not include the contribution coming from the 
displacement of the bubble in the opposite direction, i.e. for negative values of a. To include 
that one has to multiply the above expression by 2 to get the final correct expression for 

(A*))- 

(7(0) = (1(0)) = Aerfc[-^L=\. (7) 

2 V2C(0) 

Next we evaluate the correlation function (I(t)I(0)). Here also we follow the same tech- 
nique and write the correlation function as 



J^2 °9 °9 / 2ik\ ^2 X p (t)4> p (a) 2ik 2 ^ X p (0t)(f> p (a) \ 

(I(t)I(0)) = (—) J dm J da 2 J dh J dk 2 (e ^ e " =1 j 

a a > ' 

Next introducing h(s) = 2ki(f) p (a)5(t — s) + 2/c 2 0p(a)<5(s) one writes the average as 



/ 2ifci ]T X p (t)<t> p (a) 2ik 2 Y / X p (0t)(f, p (a)\ I i jP f dsX p (s)h(s)\ 

I c p=i c p=l \ — / p p=l— 00 

\ / V 

Similarly the above average can be written as 



00 00 00 



iJ2 I dsX p (s)h(s)\ -\ J dti J dt 2 Y, h (ti)H- l (ti-t2)h(t 2 ) 

, p=l-oo \ — g 



-00 — oc p— 1 



Then one gets 
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i J dsX p (s)h(s)\ 



p=l- 



-2(fc2+fc2)c( )-2fc 1 fc 2 C(t) 



where 



C(t) = ]T0 p (a) 2 (X p (t)X p (O)> (8) 
p=i 

To get the final result one has to perform integrations over fci, &2, «i and oli- Unfortu- 
nately one of the integrals over a (say a 2 ) can not be performed analytically and has to be 
carried out numerically. Now as done earlier in the calculation of (I(t)) here also one should 
consider the contribution due to the displacement of the string in the opposite direction 
which corresponds to negative a. This means not only the above expression should be mul- 
tiplied by a factor of 2 but also there would be two cross terms (9(R(a, t) — a)9(R(a, 0) + a)) 
and (9(R(a,t) + a)9(R(a,0) — a)). These cross terms would contribute equally to the cor- 
relation function. Including all these contributions, the final correct expression for the 
correlation function becomes 

</(t)/(0)> = j£ ( fll (t) + a 2 (t)) (9) 
2^2vrC(0) 

where 

. , 7 , -JL , r 2C(0)a - C(t)a 2 , . . 

aAt) = / da 2 e ^erfc — K ' w 10 

I ' 2pC{mCW-C{tf) 

and 



f , —?L , r 2C(0)a + C(t)a 2 
a 2 (t) = / da 2 e gWer/c -p 11 
-i 2 V / 2C(0)(4C(0)^-C(t)2) 

Now in principle one can use Eq. (J7|) and Eq. (JUJ) to calculate G(i) defined in Eq. (J5J). 

But to do that one has to evaluate C(t) first. Now as mentioned earlier <p p (n) introduced 

in Eq. ([6]) has to vanish at the boundaries which is satisfied by choosing 4> p (n) = sin( £ ^ ; ). 

One can show that the time dependent coefficients X p (t) obeys the following equation. 

C P ^^ = -k P X p (t) + f p (t). (12) 



8 



C P = 2L( 

and 

,2fc7rV 
k p = ( + 2kL) 

and the / p (t)'s are the random forces which satisfy 

(/,(*)> = o 

(f p {t)f p {s)) = 2( p k B T8{t-s) 
From Eq. ffT2l and using the statistical properties of the random forces one can derive 

<X p (t)X p (0)> = *j^ e -£ (13) 

where 



r, 



K p kn 2 p 2 + kL 2 



Next we will choose a = |f which corresponds to the mid point of the bubble. This further 
reduces the expression for C(t) to 



(kn 2 p 2 + KL 2 ) t 

e ^ 



C(t) = k B T £ (14) 

p=odd \ l 

Then using Eq. (113j) and replacing the sum in Eq. (TI4"1) by an integral one can write C(t) 



as 



= / *4nr^ (15) 

The factor | comes because only odd modes contribute to the sum. Fortunately the 
above integral is analytical. After carrying out the integration one gets 

C{t) = M=erfciM}. (16) 

8V KK V *5 

Using this expression for C{t) and Eq. flTJ) and Eq. (E} one can evaluate G(t) defined in 
Eq. ©. 
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III. COMPARISON WITH ALTAN-BONETT EXPERIMENT 



Here we make a comparison of our model with the experimental data obtained by Altan- 



Bonett et al 



.The fitting function they used has already been mentioned at the beginning 
( Eq. (T5])). Figure 1 is a comparison of our model with Altan-Bonett fitting function (Eq. 
02])). The comparison are made for a fixed set of parameters (a = 1, L = 10, k — 100, ( = 
5, ksT = 1) while changing the force constant of the confined harmonic well (k) which 
is also a measure of the strength of H-bonding and staking interaction between the base 
pairs. Changing k by keeping other parameters unchanged physically means changing the 
relaxation time of the normal modes of the chain as r n = r = , , \ f T , . Thus a smaller 
value of k results slow relaxation of the normal modes and a larger k results faster relaxation 
of normal modes. Another interesting observation is that the relaxation time for the higher 
normal modes (larger p) have very weak dependence on k, while the relaxation times for 
the lower normal modes (smaller p) have stronger k dependence. One can see from Fig. 1 
that there exist at least one value of k (when all the other parameters are kept fixed) for 
which a very good comparison with Altan-Bonett results can be made. Moreover a small 
change in the value of k results poor comparison with the experimental result as shown in 
Fig. 1 (dashed, blue). A smaller value of k makes the dynamics slower as expected and 
also a higher value results faster decay of the autocorrelation function (green, dashed-dot). 
This strong k dependence also suggests that the lower normal modes of the strand mostly 
contribute to bubble dynamics. Here we would also like to mention that the choice of a 
particular set of parameters is not unique as more than one set of values can also be used to 
make a reasonably well comparison as is also found with the model of Chatterjee et al jsj. 



A. long and short time behavior of G(t) 

In this section we explore the short and long time behavior of G(t). Let us first consider 
the correlation function used to fit the experimental data by Altan-Bonett Eq. (j2]). At 
short time one can approximate, erfc[z] ~ 1 — ^= as (z — > 0) and the correlation function 
becomes G ex (z) = 1 — 2^J^z 1 / 2 . Thus it behaves as a power law. To explore the long 
time limit we use the following approximation, erfc[z] ~ as (z — >■ oo), which gives 

_ az 

G ex ( z oo) = 7=7= ■ Now it would be interesting to see what happens to G(t) (or G(z), 
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z = t/ti/2) in our formulation. All time dependence in G(t) comes through (I(t)I(0)} or in 
other words through ai(t) and a 2 (t) defined earlier. To analyze the the short time behavior 
of G(t) we first rewrite short time C(t) as 

C(t) =C(0)(l-2cxy/i) (17) 

where, c x = and C(0) = 

' 1 V 8V&K 

Using the above short time expression for C(t) and keeping the leading order in t the 
short time expression of a x (t) simplifies to 

ai(t) = J da 2 e~^oJ {A 1 {a 2 ) + B x {a 2 )^/i) = K X + L x Vt (18) 

a 



with A x (a 2 ) = ^=(Vf« - 5|) and 5 i(«2) = C(0)-§(2^ Cl a 2 - f Cl a). 
Similarly 



~ a 2 

a 2 (t) = J da 2 e~^oJ(A 2 (a 2 ) + B 2 (a 2 )Vi) = K 2 + L 2 \/£ (19) 



where 

A 2 (a 2 ) = ;^=(#* + B 2 (a 2 ) = C(0)-§(-^ Cl c* 2 - f c x a), K x 



2 2 

J= 2 OO a 2 



J da 2 e 8C, (°) A x (a 2 ), L x = J da 2 e 8C (°) -Bi(a 2 ), K 2 = J da 2 e sc (°) A 2 (a 2 ), L 2 



J da 2 e sc (°) B 2 (a 2 ). 

a 

Then combining Eq. ({TBI) . Eq. ( TT9|) one finally gets the short time limit of G(t). 

G(t) = G(0)(l-A v / t) (20) 

where, A = G(0) = 1 (as G(t ) is norma li zed ) , d x = ( e r/c[^|=]) 2 , 

d 2 = . 

Hence at short time it has the similar time dependence as the one used by Altan-Bonett 
to fit the experimental data. Next we analyze the long time behavior of G(t). As mentioned 
earlier all the time dependence of G(t) comes from C(t) which is embedded in a x (t) and 
a 2 (t). As in the long time, C(t) << C(0), one can further simplifies the expressions for a x (t) 
and a 2 (t). To do this consider the complimentary error function sitting inside the integral 
in Eq. ([10]). 



11 



2C(0)a - C(t)a 2 . , r 2C(0)a - C{t)a 2 , 
evtcl — — = er/c — — — 

2y / 2C(0)(4C(0) 2 -C(t) 2 ) " 2^2C(0)(4C(0) 2 ) ' 
Now the complimentary error function is in the form erfc[P\ — Q\C(t)\ with P\ 



Qi 

~ 4 v / 2(c'(o)- 3 / 2 )' i ^ s approaches zero in the long time we can make a series expansion 
of the above complimentary error function and keep only the first two terms to get 

erfclP, - Q x C{t)) = er/c[Pi] + ^ P *g lV ^ , 

V 71 " 

Now one can analytically perform the integration over a 2 to get the long time expressions 
for ai(t). Similarly we get the long time limit of a 2 (i) by following the above steps. It shows 
that the autocorrelation function G(t) approaches zero the same way C(t) approaches zero 
in the long time limit. Thus in the long time limit G(t) approaches zero as ^ as is also 
predicted by Altan-Bonett fitting model. 

IV. CONCLUSIONS 

The stochastic dynamics of DNA bubble formed due to rupture and reformation of hydro- 
gen bonds is modeled based on a Rouse chain description of the DNA strand. Although the 
model is very simple it produces the experimental results of Altan-Bonett reasonably well. 



Unlike other well known models 



5[ our model goes beyond a single reaction coordinate or 
order parameter description by taking into account of the collective nature of the dynamics 
through different modes of the chain as is done in the Rouse description in the simplest 
possible way. The dynamics seems to be very sensitive to the relaxation times of different 
normal modes of the chain which is physically understandable as the "bubble dynamics" 
should have the contribution from all the possible normal mode of vibration of the chain. 
However the current model probably can not account for the bubble size distribution. We 
are in the process of developing a more realistic model which can account for the issue like 
bubble size distribution. 
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FIG. 1: The autocorrelation function against time (t/ii/2)- The values of other parameters used 
are: a = 1, L = 10, k = 100, £ = 5, /c B T = 1. 
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